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ABSTRACT 


We report the discovery of a massive protostar M17 MIR embedded in a hot molecular core in M17. 
The multiwavelength data obtained during 1993-2019 show significant mid-IR, (MIR) variations, which 
can be split into three stages: the decreasing phase during 1993.03—mid-2004, the quiescent phase from 
mid-2004 to mid-2010, and the rebrightening phase from mid-2010 until now. The variation of the 
22 GHz H20 maser emission, together with the MIR variation, indicates an enhanced disk accretion 
rate onto M17 MIR during the decreasing and rebrightening phases. Radiative transfer modeling of 
the spectral energy distributions of M17 MIR in the 2005 epoch (quiescent) and 2017 epoch (accretion 
outburst) constrains the basic stellar parameters of M17 MIR, which is an intermediate-mass protostar 
(M, ~ 5.4 Mo) with Mace ~ 1.1 x 1079 Mo yr-! in the 2005 epoch and Mace ~ 1.7 x 107? Mo yr7! in 
the 2017 epoch. The enhanced Moce during outburst induces the luminosity outburst AL ~ 7600 Lo. 
In the accretion outburst, a larger stellar radius is required to produce Macc consistent with the value 
estimated from the kinematics of H20 masers. M17 MIR shows two accretion outbursts (At ~ 9—20 yr) 
with outburst magnitudes of about 2 mag, separated by a 6 yr quiescent phase. The accretion ourbusrt 
occupies 8396 of the time over 26 yr. The accretion rate in outburst is variable with amplitude much 
lower than the contrast between quiescent and ourbusrt phases. The extreme youth of M17 MIR 
suggests that minor accretion bursts are frequent in the earliest stages of massive star formation. 


Keywords: Accretion, accretion disks — Stars: massive — Stars: formation — Stars: protostars — Stars: 
individual objects (M17 MIR) 


1. INTRODUCTION 


High-mass stars (M. = 8 Mo) have strong impacts on 
their environment through their various feedback affects. 
However, there is still no consensus on the basic forma- 
tion mechanism of high-mass stars. The massive star 
formation models of disk accretion and accumulated ob- 
servational evidence appear to favor the idea that high- 
mass stars form as lower-mass progenitors that continu- 
ously accrete material from their gas-rich circumstellar 
disks at high accretion rates (> 107^ Mo yr-!) (e.g., 
Hosokawa et al. 2010; Beltrán & de Wit 2016; Haem- 
merlé et al. 2016; Motte et al. 2018). Observational evi- 


dence for episodic accretion extends across a wide range 
of protostellar mass, from low values (e.g., Audard et al. 
2014, for a review) to a high mass ~ 10 — 20 Mo (e.g., 
Tapia et al. 2015; Caratti o Garatti et al. 2017; Hunter 
et al. 2017). Although the total periods of multiple 
bursting accretions constitute only a few percent of the 
entire accretion phases of massive young stellar objects 
(MYSOs), bursting accretions can contribute a substan- 
tial fraction of the zero-age main-sequence (ZAMS) mass 
of MYSOs (Meyer et al. 2019b). Moreover, bursting 
accretions influence the evolutionary tracks of MYSOs; 
during each bursting accretion, MYSOs experience rapid 
excursions toward more luminous, but colder, regions of 
the Hertzsprung-Russel diagram, which are usually oc- 
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* Based on observations collected at the European Southern Obser- 
vatory, Chile, Programs 073.C-0170(A+B), and 077.C-0174(A). 


cupied by evolved massive stars (Meyer et al. 2019a). 
This short bloating of MYSOs can occur multiple times 
up to the end of the pre-main-sequence evolution, dif- 
ferent from the spectroscopically confirmed bloating of 
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MYSOs during the contraction phase toward the main 
sequence (Chen et al. 2015; Ramírez-Tannus et al. 2017). 
The discovery of burst-induced photometric variability 
for MYSOs is impossible at optical wavelengths and dif- 
ficult in the IR, due to severe extinction and the large 
distances of MYSOs. 

At a distance of ~ 2.0 kpc (Xu et al. 2011; Chibueze 
et al. 2016; Wu et al. 2019), M17 is among the best lab- 
oratories in the Galaxy for investigating the formations 
of high-mass stars. Several efforts focused on the high- 
mass stars in M17 have attempted to measure the phys- 
ical properties of MYSOs, e.g., Hanson et al. (1997), 
Nielbock et al. (2001), Hoffmeister et al. (2006), and 
Ramírez-Tannus et al. (2017). Meanwhile detailed stud- 
ies of several peculiar MYSOs provide better constraints 
on their evolutionary stages and circumstellar environ- 
ments, e.g., Chini et al. (2000, 2004, 2006), Kassis et al. 
(2002), Nielbock et al. (2007), Chen et al. (2015), and 
Lim et al. (2020). Apart from the classified MYSOs 
(mostly IR-bright), it is very likely that unknown high- 
mass protostars are deeply embedded in dense cores 
within M17. There are nine sites of 22 GHz H2O masers 
in M17 (Johnson et al. 1998; Breen et al. 2010); however, 
only one is likely associated with the known massive pro- 
tostar IRS5A (Chen et al. 2015); the driving sources of 
the remaining H20 masers are still unclear. 

The expanding shell traced by the H20 maser spots 
(Chibueze et al. 2016, hereafter CJO16) motivates the 
search for the driving source. The coordinates of the pu- 
tative driving source are determined as R.A.(J2000) = 
18520™238017, decl.(J2000) = —16?11'47/98 with posi- 
tional errors of the order of 1 mas (CJO16), which co- 
incides with a faint point source (hereafter M17 MIR) 
with Spitzer 4.5 and 5.8 um magnitudes of 9.97 and 7.82 
mag, respectively (see more descriptions of GLIMPSE 
in Sect. 2.3). On the other hand, none of the com- 
pact radio sources reported by Rodríguez et al. (2012) 
matches M17 MIR, which eliminates the association 
with a hyper- or ultracompact H II region, and thus 
indicates a very early stage for M17 MIR. 

'This paper reports the analyses of the multiwavelenth 
data available for M17 MIR. The multiwavelength data 
are described in Sect. 2. The analyses of these data and 
the results are presented in Sect. 3. The variability of 
M17 MIR is classified in Sect. 4. The derived properties 
of M17 MIR are discussed in Sect. 5, and the conclusions 
are summarized in Sect. 6. 


2. OBSERVATIONS 


2.1. IR Imaging Data from the Very Large Telescope 


The wide-field JH K,- and L'-band images were ob- 
tained with ISAAC at the Very Large Telescope (VLT) 


in 2002 September and 2004 May to September, re- 
spectively. The pixel scale is 07148 and 0"07, and the 
FWHM is 0/6 and 0"5 for the JH K,- and L'-band im- 
ages, respectively (see more details of observations in 
Hoffmeister et al. 2008). These JH K,L'-band images 
and photometric data were published by Hoffmeister 
et al. (2008). 

'The MIR images were obtained toward the M17SW 
photodissociation region (PDR) (e.g., Crété et al. 1999; 
Pérez-Beaupuits et al. 2012) with the VLT imager and 
spectrometer for mid-infrared (VISIR; Lagage et al. 
2004) at the VLT in 2006 April through the SiC filter 
(A = 11.8 um, AA = 2.1 um). The data were taken with 
the VISIR. img obs AutoChopNod template. The photo- 
metric standard star was HD 178345 during this obser- 
vation. The data were reduced and calibrated using 
the EsoReflex workflows (Freudling et al. 2013) and the 
VISIR pipeline recipes 4.3.7! . The image is of good qual- 
ity with FWHM z 04, reaching the diffraction limit. 


2.2. The Infrared Space Observatory (ISO) Data? 


The archival imaging data of ISO (Kessler et al. 1996) 
were obtained for the M17 SW PDR from the ISO Data 
Archive?. The ISO observations for the M17 SW PDR 
were carried out in 1997 March and published in Crété 
et al. (1999); they aim at the interstellar medium in- 
side the PDR. The calibrated ISO/ISOCAM images in 
the three broad filters (LW1, LW4, and LW6) with cen- 
tral wavelengths of 4.5, 6.0, and 7.7 um are used in this 


paper. 


2.3. The Spitzer Data 


'The Galactic Legacy Infrared Mid-Plane Survey Ex- 
traordinaire (GLIMPSE; Benjamin et al. 2003; Church- 
well et al. 2009) is a legacy science program of the Spitzer 
Space Telescope that surveyed the inner Galactic plane 
in the 3.6, 4.5, 5.8 and 8.0 um bands with the IR cam- 
era IRAC. The angular resolution in the four bands 
(11,12,I3,and I4) is about 275. The GLIMPSE mosaicked 
images of M17 were obtained from the NASA/IPAC In- 
frared Science Archive^; the median observation date 
was 2005 September. The GLIMPSE point sources were 
retrieved from the GLIMPSE I Spring 07 Archive?. 


1 http://www.eso.org/sci/software/pipelines/ 


? Based on observations with ISO, an ESA project with instru- 
ments funded by ESA Member States (especially the PI countries: 
France, Germany, the Netherlands and the United Kingdom) and 
with the participation of ISAS and NASA. 

3 https:/ /www.cosmos.esa.int /web/iso/access-the-archive 


^ bttp://irsa.ipac.caltech.edu/frontpage/ 


5 See more descriptions in http://irsa.ipac.caltech.edu/Missions/spitzer.html 
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Figure 1. Multiwavelength images centered on M17 MIR. Proper motions of H2O maser spots (CJO16) are overlaid on the 


image at 4.5 um. 


M17 was revisited by Spizter in the warm phase in 
2017 July. During these observations, only the two 
bluest bands (3.6 and 4.5 um) were available, and the 
flux-calibrated images (program ID: 13117, observer: 
Robert Benjamin) were obtained from the Spitzer Her- 
itage Archive. 


2.4. The Herschel Data 


In order to cover the far-IR regime, we incorporated 
the publicly available imaging observations from the 
Herschel Space Observatory (Herschel) and its PACS in- 
strument at 70, 100, and 160 um (Poglitsch et al. 2010). 
'The Herschel data used in this work are Level 3.0 prod- 
ucts (eight observations with IDs from 1342192767 to 
1342192774) provided by the Herschel Science Archive". 
The pixel scales are 1"6, 1"6, and 32 at 70, 100, and 
160 um, and the measured resolutions are 870, 870, and 
12"0, respectively. 


2.5. The WISE/NEOWISE data 


The Wide-field Infrared Survey Explorer (WISE, 
Wright et al. 2010) is a 40 cm telescope in a low earth 
orbit that surveyed the entire sky in 2010 using four IR 
bands at 3.4, 4.6, 12, and 22 um. The angular resolution 
in the four bands (W1, W2, W3, and W4) are 6"1, 64, 
6/5, and 12"0, respectively. Only the WISE W1 and 
W2 data of M17 are considered in this work, because 
the M17 H II region is saturated in the W3 and W4 


6 Hosted at https://irsa.ipac.caltech.edu/Missions/spitzer.html 
7 http://archives.esac.esa.int/hsa/whsa/ 


bands. With the primary aim of studying near-Earth 
Objects, NEOWISE observations resumed in 2013 De- 
cember and continue to date, after the telescope's cryo- 
gen tanks were depleted. Due to the orbit of the tele- 
scope, WISE/NEOWISE observations visit M17 twice 
a year, in March and September. From the NEOWISE 
2020 data release containing W1 and W2 observations 
from 2013 December until 2019 December, the single- 
exposure data of M17 collected over 6 yr (12 epochs) 
are analyzed in this work, in addition to the WISE all- 
sky survey image of M17 (two epochs). 

The single-exposure images of the WISE/NEOWISE 
mission toward M17 were retrieved from the WISE data 
hosted on IRSA *, with the service of WISE/NEOWISE 
Coadder. The single exposures were coadded to obtain 
higher-quality W1/W2 images with enhanced resolution 
of 076875/pix. Finally WISE/NEOWISE images from 
14 epochs of M17 were used in this work. Aperture 
photometry with a fixed radius of 8.5 pixel (= 5^8) is 
applied to all WISE/NEOWISE images with a field of 
view of 3’ x 3’. Because the WISE/NEOWISE W2 band 
is very similar to the Spitzer I2 band, the photometric 
results from multi epochs are checked against the iso- 
lated star G015.0198-00.6768 located 33" south of the 
image center. G015.0198-00.6768 is not cataloged in the 
ALLWISE catalog (Cutri et al. 2021) or in the unWISE 
catalog (Schlafly et al. 2019), but has a cataloged Spitzer 
I2 magnitude of 8.55 (Spitzer Science 2009). 


8 https:/ /irsa.ipac.caltech.edu/Missions/wise.html 
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2.6. SOFIA Data 


The SOFIA? (Young et al. 2012) observations toward 
the M17 SW PDR were carried out with the FORCAST 
instrument (Herter et al. 2013) at 19 wm (Aeg = 19.7 um; 
AA = 5.5 um) and 37 um (Aeg = 37.1 um; AA = 5.5 um) 
on 2017 August 2 (PI: James M. De Buizer). The cali- 
brated Level 3 images at 19 and 37 um were downloaded 
from the SOFIA Data Cycle System !?. The resolutions 
at 19 and 37 um are 3"3 and 3"5, respectively (De Buizer 
et al. 2017). The SOFIA/FORCAST images fill the gap 
between the VLT/VISIR and Herschel images in wave- 
length coverage. 


2.7. JCMT SCUBA2 Submillimeter data 


The James Clerk Maxwell Telescope (JCMT) 
SCUBA2 submillimeter observations were taken at 
450 um and 850 um on 2016 April 10 for M17, as a 
part of the JCMT large program “SCOPE: SCUBA- 
2 Continuum Observations of Pre-protostellar Evolu- 
tion” (JCMT program code: M16AL003). The pipeline- 
reduced and calibrated maps of M17 at 450 um and 
850 wm were retrieved from Canadian Astronomy Data 
Centre !!, 


3. RESULTS 


3.1. A very early high-mass protostellar object driving 
a protostellar outflow 


The earliest detection of M17 MIR was at 3.9 um L- 
band by the Canada-France-Hawaii Telescope (CFHT) 
observations in 1992 December and 1993 May with a 
spatial resolution of 1” (Giard et al. 1994). This ob- 
ject, named as source 2 by the authors, and six other 
extremely red objects are suggested to be young stars 
in very early evolutionary stages (see Table 1 in Giard 
et al. 1994). Figure 1 presents the multiwavelength im- 
ages centered on M17 MIR. A summary of the IR flux 
densities obtained by the various facilities except for the 
earliest CFHT observation is presented in Table 1. 

The source sizes in the various IR images are in 
proportional to FWHMs of the corresponding images. 
In the 11.8um image of 0/4 resolution (= 800au), 
M17 MIR is still a single source. M17 MIR appears 
as a compact point source at 3 — 37 um. At Herschel 
70um and 1604m, the far-IR emission is dominated by 
the contamination from M17 IRS5 and M17 UCI to the 


9 SOFIA is jointly operated by the Universities Space Research 
Association, Inc. (USRA), under NASA contract NAS2-97001, 
and the Deutsches SOFIA Institut (DSI) under DLR contract 50 
OK 0901 to the University of Stuttgart. 

10 https:/ /dcs.arc.nasa.gov/ 


11 http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/jcmt/ 


'Table 1. IR flux densities 


Telescope — Aeg Epoch Flux Densities FWHM 
(um) (mJy) (arcsec) 
VLT 3.78 2004.09 1.91 + 0.17 0.5 
Spitzer 4.50 2005.07 18.4 + 2.2 2.0 
Spitzer 5.71 2005.07 85.0 + 8.0 2.0 


ISO 7.72 1997.03 (2.2 0.3) x 10? 3.7 
Spitzer 7.85 2005.07 (1.4 0.3) x 10? 2.0 
VLT 11.80 2006.04 (1.42 0.07) x 10? 0.4 
SOFIA 19.7 2017.08 (5.1 1.2) x 10? 3.3 
SOFIA 37.1 2017.08 (7.3 1.2) x 10? 3.5 
Herschel 70 — 2010.03 < 7.7% 10 8.0 


NoTE—Central wavelengths and effective bandwidths are 
adopted from the SVO Filter Profile Service (Rodrigo & 
Solano 2020). 


northeast. M17 MIR is not resolved at its IR location 
of M17 MIR, likely due to the strong ambient far-IR 
emission. For M17 MIR, we estimated an upper limit 
< 77 Jy at 70 um. 

The faint IR brightness of M17 MIR suggests that it is 
a low-luminosity object or the extinction to M17 MIR is 
quite high. From a quick examination of the IR images 
of M17 MIR, we think high extinction is the more likely 
reason. M17 MIR’s brightness shows a steep rise in the 
3 — 8 um range (see also Figure 8). However, M17 MIR 
is not visible in the ISO broad filter image at 9.6 jum 
(see Figure 2a-f; Crété et al. 1999). We know that the 
silicate absorption feature at 9.7 jum can cause high and 
broad absorption to a highly obscured object at around 
this wavelength. The high extinction to M17 MIR can 
explain the nondetection at 9.6 um. 

The IR flux densities of M17 MIR are at least two 
orders of magnitude lower than those of the IR-bright 
high-mass protostar IRS5A (Chen et al. 2015) within 
the same M17 SW PDR. On the other hand, the out- 
wards motion of the CJO16 H20 masers very likely 
traces the shocked gas accompanying the protostellar 
outflow driven by M17 MIR. In the earliest stage of 
high-mass protostars, the mass loss might evolve from 
a highly uncollimated expanding shell to a collimated 
outflow (e.g., Carrasco-González et al. 2015). The dy- 
namical age of the expanding shell driven by M17 MIR 
is 12.5 yr (CJO16), younger than the dynamical age 
(~ 25 yr) of W75N(B)-VLA 2 (Carrasco-González et al. 
2015). The faint IR brightness of M17 MIR and the 
young dynamical age of the CJO16 H20 masers suggest 
that M17 MIR is a promising candidate for a high-mass 
protostar in the earliest stage. The analysis of the spec- 
tral energy distribution (SED) constructed from Table 1 
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Figure 2. JCMT-SCUBA2 850 um map of M17 SW in left, and the HCO+ J = 9 — 8 integrated intensity map of M17 SW in 
right (Pérez-Beaupuits et al. 2015) . The 22 GHz water maser sources detected by Johnson et al. (1998) are denoted by crosses 
with numbers as same as Johnson et al. (1998). The white circle in both maps denote the position of M17 UCI. 


will help to better constrain the fundamental properties 
of M17 MIR (Sect. 5.3). 

'The warm area close to the central object emits radia- 
tion mostly at MIR wavelengths, while the surrounding 
envelope of much larger size is expected to radiate at 
submillimeter to millimeter wavelengths due to a much 
lower effective temperature. Figure 2 presents maps 
of the 850m dust continuum emission and the HCO+ 
J —9— 8 line integrated intensity of M17 SW PDR. 
M17 MIR, represented by the H20 maser source 2 de- 
tected by VLA observations (Johnson et al. 1998, here- 
after JDG98), is located exactly in the secondary peak 
of the 850 um dust continuum emission. This dust con- 
tinuum emission peak, along with the other two peaks, 
make up the northern condensation in M17 SW PDR 
that was first noticed by Stutzki & Guesten (1990). 
The three main components in the northern condensa- 
tion in M17 SW PDR were studied in detail through 
submillimeter and millimeter observations. Their phys- 
ical properties were estimated by Hobson et al. (1993) 
(their Table 4). M17 MIR is located within the FIR3 
source in Hobson et al. (1993). On the other hand, the 
map of integrated line intensity (range: 0 — 40km s7}, 
peak velocity ~ --20kms^!) of the HCO+ J = 9 > 8 
line shows clumpy structures similar to the dust con- 
tinuum emission. A dense core whose size is compara- 
ble to the beam size of the observations (c 8", or 0.08 
pc) is coinciding with M17 MIR. Given the critical den- 
sity nei, = 9 x 107 cm7? and an excitation temperature 
Eup = 192.58 K of the 802.458 GHz HCO+ J = 9 > 8 
line (Pérez-Beaupuits et al. 2015), the dense core associ- 


Table 2. L-band magnitudes of embedded YSOs in M17 SW 
during the 1993.03 and 2004.09 epochs. 


Name R.A. Decl. 1993.03 2004.09 
(J2000) (J2000) (mag) (mag) 
2 = M17 MIR 18:20:23.02 -16:11:47.9 10.2 12.78 
3 18:20:23.91 -16:11:43.4 10.8 10.38 
5 18:20:24.10 -16:11:40.1 11.2 10.99 
T 18:20:22.34 -16:11:28.7 8.8 8.76 


NoTE-—See Figure 3 for the positions of these sources 
on the Spitzer I1 and I2 images. 


ated with M17 MIR is likely a candidate for a hot molec- 
ular core (HMC; Cesaroni 2005). HMCs commonly cor- 
responds to a hydrogen column density Ny = 10?? or 
even higher (e.g. Bonfand et al. 2017). Although a high- 
mass protostar embedded in an HMC is undetectable at 
short wavelengths due to high extinction, the detection 
of the central protostar becomes possible at MIR wave- 
lengths. For instance, the IR counterpart F4 of the HMC 
G9.62+0.19-F appears at 3.8 jum in the L’ band and at 
longer wavelengths (Linz et al. 2005). 

'The multiwavelength data of M17 MIR and the sub- 
millimeter to millimeter data of its surrounding medium 
indicate that M17 MIR is a very early high-mass proto- 
star embedded in an HMC and is also driving a proto- 
stellar outflow. 


4. VARIABILITY 
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4.1. Infrared variability 


For M17 MIR, the photometric magnitude in the 1993 
CFHT observation was 10.2 + 0.2 (source 2 in Table 2; 
Giard et al. 1994), while the magnitude in the 2004 VLT 
observation was 12.78 + 0.14 (Hoffmeister et al. 2008). 
The filter profiles in the two observations are almost 
identical. Table 2 lists the L-band magnitudes of the 
embedded YSOs in the two observations. Interestingly, 
M17 MIR is about 10 times fainter in the 2004.09 epoch 
than in the 1993.03 epoch, whereas the other sources are 
fairly stable. 

The VLT L’-band observations in 2004.09 were one 
year earlier than the Spitzer observations in 2005.09. By 
interpolating the Spitzer flux densities at I2 and I3 to 
that at L’ under the assumption of a blackbody radia- 
tion, the flux density of M17 MIR in the L’ band is about 
2.2mJy in the 2005.09 epoch. This predicted flux den- 
sity in 2005.09 is close to the 1.91 + 0.17 mJy measured 
in 2004.09. We therefore speculate that M17 MIR's lu- 
minosity is stable during 2004.09-2005.09. 

M17 MIR shows variabilities in the Spitzer I1 and I2 
images from 2005.09 and 2017.07 (see Figure 3). Com- 
paring the flux densities at I1 in the two epochs for an 
aperture radius of 16, we find an increase by a factor of 
2.9. Similarly, the flux densities at I2 show an increase 
by a factor of 3.3 in an aperture radius of 3/2. In con- 
trast, the other sources are stable in the two bands (see 
Table 3). 

'The MIR variabilities of M17 MIR found for different 
epochs motivate us to compare its flux densities from 
the ISO observations in 1997.03 and the Spitzer obser- 
vations in 2005.09 and 2017.07. Table 4 shows the flux 
densities of M17 MIR and other nonvariable sources in 
the ISO/LW1 and Spitzer/I2 bands, which have almost 
identical filter profiles centered at 4.5 um. Figure 4 com- 
pares the ISO LW1 and Spitzer I2 images obtained in 
1997.03 and 2005.09, respectively. Direct comparison 
from Figure 4 clearly indicates that M17 MIR is signifi- 
cantly brighter than sources 3 and 5 in the ISO LW1 im- 
age in the 1997.03 epoch, while it becomes fainter than 
sources 3 and 5 in the Spitzer I2 image in the 2005.09 
epoch. Given that sources 3 and 5 are unchanged (see 
context above), M17 MIR is varying between 1997.03 
and 2005.09. The aperture photometry of M17 MIR in 
the 4.5 um band shows variability from 1997 to 2017. 

For M17, WISE/NEOWTSE images in 14 epochs, in- 
cluding 2010.03 and 2010.09, and March and September 
in 2014 to 2019, were analyzed in this work. M17 MIR is 
only detected in the W2 band. The WISE W1 flux den- 
sities, extrapolated from its flux densities in the Spitzer 
Il and I2 bands, are about 0.6 and 1.8 mJy during the 
2005.09 and 2017.07 epochs, respectively. The extreme 
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Figure 3. Top panel: Spitzer/Il images in 2005.09 (left) 
and 2017.07 (right). Bottom panel: Spitzer/I2 images in 
2005.09 (left) and 2017.07 (right). The embedded YSOs 3, 
5, and 7 are the same as in Table 2. 
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Figure 4. ISO/LW1 image in 1997.03 (left) and Spitzer/I2 
image in 2005.09 (right). The embedded YSOs 3, 5, and 7 
are identical to those in Figure 3. 


faintness at W1 is consistent with the nondetection of 
M17 MIR in the corresponding W1 images. 

Figure 5 presents the WISE/NEOWISE W2 images 
centered on M17 MIR in 14 epochs. The pixel resolution 
of 076875 clearly separates M17 MIR from the nearby 
IR source Anon 3 to the southwest. The photometric re- 
sults of M17 MIR, source 7, and G015.0198-00.6768 are 
tabulated in Table 5 for 14 epochs. In 2010.03, the W2 
flux density of M17 MIR was 17.7 + 4.4 mJy, which is 
very close to the value 18.4+2.0 in the Spitzer I2 band in 
2005.09. Given the fact that the WISE/NEOWISE W2 
band is very similar to the Spitzer I2 band, M17 MIR 
was roughly stable in brightness from 2004 to 2010. 
However, six months later M17 MIR was twice as bright 
in 2010.09, with a flux density of 33.7 + 4.4 mJy. It is 
very likely that the rebrightening of M17 MIR started 
between 2010.03 and 2010.09. In the observations of 
NEOWISE, M17 MIR generally shows a steady rise in 
brightness, but with a jump and fall motion epoch by 
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Figure 5. WISE/NEWOWISE W2 images centered on M17 MIR in 14 epochs. 


epoch in its light curve (see Figure 6). Between 2016 
and 2017, M17 MIR underwent a steep rise in brightness 
with an amplitude of 1.7. In 2018 M17 MIR’s brightness 
dropped back to the level of 2016. and was followed by 
a rise with a mean amplitude of 1.5 between 2018 and 
2019. 


4.2. Combined light curve in the L/ and W2 bands 


'The MIR flux measurements of M17 MIR span several 
wavelengths and observation epochs. The 3.8 um L’ fil- 
ter is closest to the L filter of the 1993.03 epoch observa- 
tion (Giard et al. 1994), while the observations of other 
epochs provided flux densities in at least two filters. For 
a direct comparison from 1993 to 2017, the flux density 
in the 3.8 um L’ filter is estimated for 1993.03, 1997.03, 
2005.09, and 2017.07. The multiwavelength flux densi- 
ties of M17 MIR in 1997.03, 2005.09, and 2017.07 epochs 


yield a reasonable estimate for the continuum in the 
range 3 — 8 um, assuming blackbody emission. The flux 
density of M17 MIR in the 3.8 jum L’ filter is the con- 
tinuum emission inside the transmission curve of the L’ 
filter. Note that the blackbody temperatures for the 
three different epochs are different. For the L-band ob- 
servations in 1993, the three blackbody temperatures 
are used for estimating the L’ flux density of M17 MIR. 
We find that the largest difference among the three val- 
ues is less than 1596, thus the mean value is used as the 
L/ filter flux density. For 1993.03, 1997.03, 2005.09, and 
2017.07 we obtain L’ flux densities; for 2004.09 there 
is direct L/-band photometry for M17 MIR. Figure 6 
presents the combined light curve of M17 MIR in the L’ 
and W2 bands. The time span of the multiepoch data 
is 26.5 yr. 


8 CHEN, ZHIWEI ET AL. 


Table 3. Flux densities of M17 MIR Measured from the 
Spitzer Observations in 2005 and 2017 


Name Spitzer 2005 2017 Aperature 
Band (mJy) (mJy) (arcsec) 
M17 MIR I1 1.0 + 0.4 2.9 + 0.5 1.6 
I2 19.0 + 2.5 52.9 + 4.7 3.6 
3 I 9.7 X: 0.9 9.6 + 0.8 2.0 
I2 33.8 + 3.7 33.7 + 3.9 2.5 
5 Il 10.7 +0.7 8.3 + 0.6 2.0 
I2 23.8 + 1.5 20.9 + 1.3 2.0 
7 I1 37.8 + 5.0 31.83: 3.7 2. 
I2 201.2 +17.2 193.0 + 6.9 4.0 


NoTE—Aperture photometry on the sources in this 
table uses fixed in positions, aperture size, and circular 
annulus size in the 2005 and 2017 epochs, for direct 
comparison between the two epochs. 


Table 4. Flux densities of M17 MIR at 4.5 um measured 
with ISO/ISOCAM in 1997 and Spitzer in 2005 and 2017. 


Name 1997 2005 2017 Aper. 
(mJy) (mJy) (mJy) “) 
M17 MIR 32.1 11.1 19.0 + 2.5 52.9 + 4.7 3.6 
T 202.7+ 31.9 201.2 +17.2 193.0 + 6.9 4.0 


NoTE—The aperture photometry of M17 MIR in the 
ISO/LW1 image might be affected by the row of bad 
pixels lying below M 17 MIR, thus its flux density 
might be underestimated. 


M17 MIR shows decreasing, quiescent and rebrighten- 
ing phases in its MIR light curve, as shown in Figure 6. 
Due to the lack of 3— 4 um observations before 1993 and 
during 1997—2004, the duration of the decreasing phase 
is poorly constrained. From the available data, the qui- 
escent phase between the decreasing and rebrightening 
phase might last about 6 yr, from mid-2004 to mid-2010. 
Because M17 MIR seems to continue the trend of re- 
brightening, the timescale of the rebrightening phase is 
Z 9 yr. The continuing NEOWISE observations will 
provide vital information on the rebrightening timescale. 


4.3. H20 Maser Variability 


Table 6 summarizes the 22 GHz H2O maser observa- 
tions from 1984 to 2008. The earliest H3O maser de- 
tection associated with M17 MIR was achieved at the 
VLA in C configuration in 1984 June (Forster & Caswell 
1989, hereafter FC89). The complete catalog of H20 
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Figure 6. Top: light curve of M17 MIR at 3.8 um in the L’ 
band; Bottom: light curve at 4.6 um in the W2 band. The 
time span is over 26 yr from 1993 to 2019. 


Table 5. WISE/NEOWISE flux densities of M17 MIR and 
other IR sources 


Epoch MJD M17 MIR source 7 G015.0198 @ 
W2 (mJy) W2 (mJy) W2 (mJy) 
2010.03 55,280.5 17.6 (4.4) 250.0 (5.1) 70.7 (2.4) 
2010.09 55,463.5 33.7 (4.4) 257.8 (5.7) 72.4 (2.3) 
2014.03 56,745.5 50.2 (4.1) 210.3 (6.4) 71.9 (2.3) 
2014.09 56,9246 82.0 (4.1) 180.5 (6.4) 72.0 (2.4) 
2015.03 57,107.5 62.0 (4.4) 177.3 (6.4) 69.9 (2.4) 
2015.09 57,284.5 81.7 (4.6) 203.3 (7.7) 69.0 (3.2) 
2016.03 57,471.5 71.5 (4.7) 180.2 (7.5) 68.9 (2.8) 
2016.09 57,645.5 82.5 (4.5) 212.8 (7.7) 66.3 (3.3) 
2017.03 57,838.5 95.6 (5.3) 199.4 (7.7) 69.2 (2.9) 
2017.09 58,006.5 126.0 (4.2) 193.4 (6.3) 72.2 (2.4) 
2018.03 58,202.6 72.7 (4.3) 195.0 (5.8) 67.3 (2.4) 
2018.09 58,366.5 110.5 (4.4) 189.2 (6.3) 70.0 (3.5) 
2019.03 58,566.5 118.0 (3.9) 185.7 (5.1) 74.9 (1.6) 
2019.09 58,730.5 95.7 (4.2) 228.5 (7.8) 68.3 (2.4) 


"The Spitzer I2 magnitude of this source is 8.55, 
corresponding to 68.05 mJy, while its 
WISE/NEOWISE W2 mean flux density is 70.49 mJy. 


masers found by Forster & Caswell (1989) was compiled 
by Forster & Caswell (1999). Three FC89 H20 masers 
are associated with M17 MIR, and show two velocity 
components of 16 — 18kms~! and 23.3kms-!. One 
H20 maser coinciding with M17 MIR was detected with 
the VLA in D configuration and shows a flux density 
of 18.57 Jy (Massi et al. 1988, hereafter MCF88). One 


BURSTING MASSIVE PROTOSTAR M17 MIR 9 


'Table 6. Water masers observations 


Epoch Obs. R.A. Decl. Sep. V peak S peak Beam Ref. 
(J2000) (J2000) 


(1) (2) (3) (4) (5) (6) (7) (8) (9) 


(arcsec)  (kms-1) (Jy) 


1984.06 VLA-C 18:20:22.91  -16:11:48.29 Lj 23.32 1.5 471 x 272 FC89 
1984.06 VLA-C 18:20:23.12  -16:11:48.29 1.5 18.05 6.7 4/1 x 272 FC89 
1984.06 VLA-C 18:20:23.05 — -16:11:48.20 0.5 16.73 6.9 471 x 272 FC89 
1984.08 VLA-D 18:20:23.07 -16:11:49.2 17.4 18.57 5/6 x 7/1 MCFS88 
1996.09 VLA-D 18:20:22.98 — -16:11:47.97 0.6 15.8 121 5/0 x 3/4  JDG98 
1996.09 VLA-D 18:20:22.98 -16:11:47.97 0.6 27.7 17 5/0 x 3⁄4 JDG98 
2003.10 ATCA-EW352  18:20:23.04 -16:11:48.4 0.6 19 197 8" x 29" B10 
2003.10 ATCA-EW352  18:20:23.04 -16:11:48.4 0.6 24 60 8" x 29" B10 
2004.07 ATCA-H168 18:20:23.04 -16:11:48.4 0.6 20 67 13" x 9" B10 
2004.07 | ATCA-H168 18:20:23.04 -16:11:48.4 0.6 24 17 13" x 9" B10 
2008.08 ATCA-6B 18:20:23.02 -16:11:47.8 0.2 21 4.2 1.7 x 0.5 BE11 


JDG98 H30 maser (source 2 in JDG98) with two veloc- 
ity components is detected around the same source. The 
velocities of the FC89 H20 masers are consistent with 
those of the two components of the JDG98 H20 maser 
source 2. Both the FC89 and JDG98 H20 masers are 
likely pumped by the interaction between the outflow 
from M17 MIR and the surrounding gas. Interestingly, 
the flux density of the JDG98 H2O maser source 2 in 
the 1996.09 epoch is ten times higher than that of the 
FC89 H2O maser in 1984.06. 

Several years after the last VLA observations, three 
water maser observations were taken by ATCA during 
three epochs. A fast fading of the water maser inten- 
sity from 197 Jy in 2003.10 (Breen et al. 2010, here- 
after B10), through 67Jy in 2004.07 (B10), to 4.2Jy 
in 2008.08 (Breen & Ellingsen 2011, hereafter BE11) 
was observed. In the earlier two epochs, the spectra of 
water maser emission show two components- a stronger 
blueshifted one at ~ 19 — 20kms~! and a weaker red- 
shifted one at 24km s^ !, which are consistent with the 
VLA observations in the earlier epochs. In the 2008.08 
epoch, two components were observed in the spectrum — 
a weaker blueshifted one at ~ 18km s^! and a stronger 
redshifted one at 21 km s71. 

'The latest water maser observations were conducted 
by very long baseline interferometry (VLBI) with VLBI 
Exploration of Radio Astrometry (VERA) in 2009.10— 
2010.12 and consist of six epochs (CJO16). The VERA 
observations detected 203 water maser spots in these 
observations. The dominant blueshifted maser spots 
(Vis, £ 22kms~*) are mostly distributed to the south 
of M17 MIR, while a few redshifted maser spots (Vis, Z 
22kms~') are distributed north to the north of it (see 


Figure 7). This feature is consistent with the stronger 
blue-shifted component and the weaker red-shifted com- 
ponent observed by the VLA and ATCA. 


5. DISCUSSION 


5.1. The Connection between the IR Variability and 
H2O Maser Emission Variability 


The IR variability of M17 MIR overlaps with the sur- 
rounding H20 maser variability in time. The connection 
between them helps to clarify the physical origin of these 
variabilities. The flux density of the JDG98 H20 maser 
source 2 in 1996.09 is ten times higher than that of the 
FC89 H2O maser source in 1984.06. This H2O maser 
variation from 1984 to 1996 and the decreasing phase 
in the IR light curve during 1993 and 2004 suggest that 
the brightness peak likely occurred between 1984 and 
1993. Moreover, the fast fading of the H20 maser emis- 
sion from 2003 to 2008 and the quiescent phase in the IR. 
light curve during 2004 and 2010 together constrain the 
occurrence of the quiescent phase to between mid-2004 
and mid-2010, characterized by both low IR brightness 
and H20 maser emission. 

VERA observations in the earliest two epochs 
(2009.10 and 2010.01) were made in the quiescent phase, 
while the later observations in the epochs 2010.04, 
2010.09, 2010.11, and 2010.12 took place in the re- 
brightening phase. The proper motions of the H20 
masers spots are traced in more than three epochs (J. 
O. Chibueze 2021, private communication). The proper 
motions of the H2O maser spots are measured in the 
transient phase between quiescent and rebrightening 
phases from 2009.10 to 2010.12. The 3D motions of the 
H20 maser spots detected by the VERA observations 
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show an expanding bubble originating from M17 MIR. 
Together with the rebrightening, we suggest that both 
phenomena are due to the accretion/outflow activity of 
M17 MIR. 

Similar flaring of the H20 masers associated with 
the outflow from the high-mass protostar NGC 6334I- 
MM1B shows a mean factor of 6.5 before and after the 
burst of MM1B in the submillimeter range (Brogan et al. 
2018). Those authors attribute the enhanced accretion 
rate to the coincident flaring in the submillimeter contin- 
uum and the H20 maser emission of NGC 63341- MM1B. 
For M17 MIR, we suggest an accretion outburst to ex- 
plain both the MIR variations and the correlated H20 
maser variations. Accompanied by the expanding mo- 
tions of H20 maser spots measured during the tran- 
sient phase from quiescent to rebrightening and the close 
positive relation between mass loss and accretion rate, 
the most promising explanation for the outburst in MIR 
brightness since 2010.03 is an enhanced accretion rate. 


5.2. Outflow rate and disk accretion rate of M17 MIR 


The proper motions and radial velocities of the CJO16 
H20 masers are overlaid on the VLT L’ image, as shown 
in Figure 7. The 3D motions of the H20 maser spots 
can be used to estimate the momentum rate of the gas 
flow, given the assumption that all the momentum of the 
ejected gas is transferred to the surrounding molecular 
environment. The momentum rate P is then estimated 
from the relation (Goddi et al. 2011) 


P = 1.5 x 107? Và, R25, (2/40) ng Mo yr kms“, 


where Vio is the mean maser velocity in units of 
10kms~!, Rioo is the average distance of H5O masers 
from the driving source in units of 100 au, 2 is the solid 
angle of the gas flow, and ng is the gas volume den- 
sity in units of 105 cm-?. Rigo is calculated as 0.85 
by averaging the distances of the H20 masers from 
the dynamical center, which is determined at the po- 
sition of R.A.(J2000) = 18"20™238017, decl.(J2000) = 
—16°11'4798 (CJO16). The gas density is about 
105 cem^? according to the discussion in Sect. 3.1. Q is 
about 4r for an expanding bubble. The above phys- 
ical properties from observations yield a momentum 
rate P ~ 4 x 107? Mo yr-! kms-! for the protostel- 
lar outflow of M17 MIR. This is comparable to out- 
flows of 107? — 10° Mo yr! kms"! driven by MYSOs 
(Moscadelli et al. 2016), and one magnitude higher than 
the outflow from a low-mass YSO in IRAS 20231--3440 
(Ogbodo et al. 2017). 

'The protostellar outflow provides an indirect measure 
of the disk accretion rate. The outflow mass rate is de- 
rived by dividing the outflow momentum rate P by the 
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Figure 7. The proper motions and Visr of the CJO16 H20 
maser. The background grayscale image is the VLT L’ im- 
age. 

outflow velocity. Given P ~ 4 x 1073 Mọ yr7!kms7! 
estimated in this work and Vio © 1.9kms-! deter- 
mined by CJO16, the outflow mass rate Mout is ~ 
2 x 107^ Mo yr-!. Assuming a ratio of outflow mass 
rate to disk accretion rate of 30% (Beltran & de Wit 
2016), the disk accretion rate onto M17 MIR is Mace ~ 
7 x 1074 Mo yr™t. 

This disk accretion rate is estimated from the H20 
maser observations during 2009 to 2010, when M17 MIR 
just started to rebrighten. The MIR emission between 
2010.09 and 2010.03 varies by a factor of about 2 in the 
WISE W2 band. Johnstone et al. (2013) simulate the 
MIR flux variation induced by an enhanced disk accre- 
tion rate and find that the MIR flux variation is propor- 
tional to the disk accretion rate. The W2 flux between 
2010 and 2019 varies by an average factor of 3. We spec- 
ulate that the disk accretion rate in 2019 is likely 3 times 
the rate estimated above, reaching ~ 2 x 107? Mọ yr“. 

5.3. Modeling of the SED in Quiescent and 
Rebrightening phases 


During the quiescent phase from mid-2004 to mid- 
2010 (hereafter the 2005 epoch), we collect IR flux 
densities at six wavelengths to construct the SED of 
M17 MIR. During the rebrightening phase, the flux den- 
sities obtained in 2017.07 (Spitzer I1 and I2 bands) and 
2017.08 (SOFIA/FORCAST 19 um and 37 um bands) 
are separated by 14 days. The NEOWISE W2 flux den- 
sity was obtained at the end of 2017 September, nearly 
two months later than the Spitzer and SOFIA observa- 
tions. We collect only the Spitzer and SOFIA observa- 
tions in 2017.07/08 to construct the SED of M17 MIR 
during the outburst phase (hereafter the 2017 epoch). 
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Figure 8. SEDs of the R17 sp-s-i models that satisfy x? — XŽest < 3ndata for M17 MIR. 


The flux densities of M17 MIR in the 2005 and 2017 
epochs are tabulated in Table 7. The 30 upper limit of 
the Herschel 70um image is 77 Jy for M17 MIR, which 
is considered as the upper limit in both the 2005 and 
2017 epochs. The additional input parameters besides 
the flux densities are distance and foreground extinction 
Ay. The distance to M17 MIR is fixed to 2.0kpc. We 
assume a wide range of 10-200 mag for the foreground 
extinction Ay. 

We employed the SED model grid for YSOs that spans 
a wide range of evolutionary stages, from the youngest 
deeply embedded protostars to pre-main-sequence stars 
with debris or no disk, and which covers a wide and uni- 
form region of parameter space (Robitaille 2017a, here- 
after R17). We fit the SED of M17 MIR using several 
R17 model sets of different model components, using the 
Python-based fitting tool (1.3 version; Robitaille 2017b). 
Based on the prior knowledge of M17 MIR, it is a highly 
embedded massive protostar with an accretion disk. The 
model sets we used are sp-s-i,sp-smi,spu-smi,spubsmi. 
All four model sets have a passive disk with an inner 
radius equal to the dust sublimation radius. sp-s-i is the 
simplest model set including only the central star and 
passive disk. The other three model sets include ambi- 
ent gas (sp-smi), envelope+ambient gas (spu-smi), and 
envelope+ambient+cavity (spubsmi). 

Fitting the SED of M17 MIR to the individual mod- 
els of the four model sets returns a x? value for each 
individual model. The criterion x? — Mn < 3Ndata Can 
generally split good and bad fits. We used this criterion 
to select good fits from the four model sets. Among 
these good fits, we use the following criteria to pick rea- 
sonable models. 

Disk mass > 0.01Mọ. The disk accretion of M17 MIR 
estimated from H20 maser motions is 7 x 1074 Mọ yr-!. 


To sustain such a high accretion rate, we require this 
criterion. 

Disk inclination angle > 30°. Because the blue- and 
redshifted H20 maser spots are distinctly separated in 
the vicinity of M17 MIR, the disk inclination angle 
should not be small. 

Exclude unphysical models. Some models are unphys- 
ical with too small radii but too high temperatures. We 
compute the Stefan-Boltzmann luminosity (4x R20 T2) 
from the radius and temperature of each model. These 
unphysical models are below the ZAMS track in the H- 
R diagram. We keep only the models that are lying 
on/above the ZAMS track in the H-R diagram. 

'The above criteria for the four model sets return rea- 
sonably good fits for the SEDs of M17 MIR in the 
2005 and 2017 epochs. Because of the degeneracy of 
the four model sets, one has to score these sets not 
only based on x? values. In order to assess the model 
sets in comparison with each other in a statistically 
robust way, Robitaille (2017a) suggested calculating 
P(D|M) œ Ngooa/N, where Ngooa is the total number 
of models from a given model set, and N is the total 
number of models in that set. Towner et al. (2019) used 
the R17 model sets to find the best YSO models for 
12 MYSOs based on the overall assessment of x? and 
P(D|M) scores. For the good fits of M17 MIR, we show 
the x2,4, and P(D|M) scores of the four model sets in 
Table 8. In the 2005 and 2017 epochs, sp-s-i is the ” best 
representation” model set according to both the x2.,, 
and P(D|M) scores. 

Among the sp-s-i model set, 44 and 13 models remain 
for the SEDs of M17 MIR in the 2005 and 2017 epochs, 
respectively. In the 2017 epoch, the SOFIA observa- 
tions at 19.7um and 37.1um are very useful in distin- 
guishing the various models of a given model set (Fig- 
ure 8). Among the 13 models in the 2017 epoch (Ta- 
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Table 7. Flux densities of M17 MIR in the 2005 and 2017 epoch 


Epoch Spitzer VLT Spitzer Spitzer Spitzer VLT? SOFIA? SOFIA? Herschel 
State 3.6um 3.78um 4.5um 5.71pm 7.8m 11.84m 19.71m 3T.lum 70pm 
(mJy) (mJy) (mJy) (Jy) (Jy) (Jy) (Jy) (Jy) 
2005 2.90.21 2132 2.9 85.0480 0.140.038 0.13940.007 «0.515012 «73-12 «77 
2017 3.8 + 0.7 59.7+5.9 > 0.139 + 0.007 0.51 + 0.12 7.3 X 1.2 < 77 


“The flux density at 11.8um in the 2005 epoch is considered as the lower limit in the 2017 epoch. 


b'The flux densities at 19.7um and 37.1um in the 2017 epoch are considered as the upper limits in the 2005 epoch. 


©The measurement at 70um in the 2010.03 epoch is considered as the upper limit. 


Table 8. x? and P(D|M) scores of R17 model sets 


2005 epoch 2017 epoch 
Model x? P(D|M) Model x? P(D|M) 
Sp-s-i 9.034 0.0044 Sp-s-i 2.41 0.0013 
sp-smi 8.241 0.0022 sp-smi 15.99 0.0017 
spu-smi 14.54 0.0012 spu-smi 14.46 0.0009 
spubsmi 9.42 0.00065 — spubsmi 11.14 0.0005 


ble 9), two models correctly reproduce the SED slope 
between 19.7um and 37.1jm. Indeed, the x., model 
is different from the second fit only in the disk inclina- 
tion angle. Aided by the observations at 19.7um and 
37.1um in the 2017 epoch, the x2... model is distinctly 
better than the other sp-s-i models. 

In the 2005 epoch, the 44 sp-s-i models show a degen- 
eracy in the long-wavelength SED beyond 124m, due 
to the lack of observations at long wavelengths. Povich 
et al. (2011) found similar SED degeneracy at long wave- 
lengths for a deeply embedded YSO with hard X-ray 
emission (PCYC source 699 from Povich et al. (2011)) 
in the Carina nebula region, when they fitted the SED 
of this object to YSO model SEDs from Robitaille et al. 
(2006). In Figure 9 of Povich et al. (2011), the distri- 
butions of the well-fit YSO models show two distinct 
groups: intermediate-mass YSOs with relatively high 
circumstellar extinction and low-mass YSOs with lower 
circumstellar extinction. The group with high circum- 
stellar extinction is consistent with the high absorption 
of the X-ray source (Povich et al. 2011). For the 44 
sp-s-i models of M17 MIR in the 2005 epoch, two dis- 
tinct families are found in the plot of infrared lumi- 
nosity integrated between 12 and 500um versus disk 
maximum radius Rsk; high infrared luminosity with 
Risk > 1000au and lower infrared luminosity with 
Rsk < 1000au. The R25* > 1000au family is more 
consistent with the general understanding of MYSOs, 
which are luminous at long infrared wavelengths and 
have large circumstellar disks. Moreover, the XE asi sp- 


-s-i model in the 2017 epoch also suggests a large cir- 
cumstellar disk for M17 MIR. The R@sk > 1000au 
family of the sp-s-i models is preferred for the SED of 
M17 MIR in the 2005 epoch. Among the nine sp-s- 
i models with R25* > 1000au in Table 9, two mod- 
els with Teff = 24,540K likely produce significant 22 
GHz continuum emission that could be detected by 
the ATCA observations during 2003-2008 (B10,BE11). 
However, no 22 GHz continuum source was reported 
around M17 MIR according to the ATCA observations. 
We excluded the two very hot sp-s-i models in the fol- 
lowing analysis. On the other hand, stellar luminosity 
in the 2005 epoch should be much lower than that in the 
2017 epoch. This constraint leads to a final collection 
of only three reasonable sp-s-i models. For the three 
sp-s-i models, we calculated the mean (Tem), (Ra), (Lx), 
and foreground (Ay), weighted by yz. These values are 
listed in Table 10 for the 2005 epoch. 

In the above analysis for the R17 models satisfying 
x? — Xl est < 3Ndata; no constraint on the foreground 
extinction Ay was applied. The good SED models in 
the 2005 and 2017 epochs converge on a narrow range 
of foreground extinction Ay — 120 — 150 mag, despite a 
broad Ay range of 10— 200 mag used in the SED fitting 
procedure. This narrow Ay range still holds for the 
other “bad” SED models. The observations at 7.8m 
and 11.8um in the 2005 epoch that cover the broad ab- 
sorption feature at 9.7um yield a strong constrain on the 
foreground extinction. This is consistent with the deep 
absorption feature at 9.7 um caused by the large ab- 
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sorption of silicate grains, which is commonly observed 
for deeply embedded high-mass protostars (e.g., Morales 
et al. 2009; Pitts et al. 2018). 


The difference in luminosity between the 2005 and 
2017 epochs indicates the luminosity burst of AL — 
7600+ 900L5 during an accretion outburst of M17 MIR. 
If AL is entirely attributed to increased accretion, we 
can predict the flux density of M17 MIR in the Spitzer 
Il and I2 bands during the outburst. We assume that 
the blackbody temperature of AL during outburst is 
1000 K. Reddenning the blackbody radiation from AL 
with a certain amount of extinction essentially matches 
the flux density in the Spitzer I1 and I2 bands. We as- 
sume a total extinction to M17 MIR of 150 mag, which 
is the sum of a foreground extinction of 130 mag accord- 
ing to the good SED models for M17 MIR and an addi- 
tional extinction of 20 mag from the circumstellar disk. 
The blackbody radiation from AL = 7600 4 900 Lo suf- 
fers from a total extinction of 150 mag, equivalent to 
0.066 x 150 mag and 0.052 x 150 mag in the Spitzer I1 
and I2 bands. This toy model yields flux density of 5.8 
mJy and 52.5 mJy in the Spitzer I1 and I2 bands, re- 
spectively, which are comparable to the observed values 
3.8 + 0.7 mJy and 59.7 + 5.9 mJy in the 2017 epoch, 
differing by a factor of less than 2. 

The luminosities of M17 MIR in Table 10 are indeed 
the total luminosites, i.e., the combination of photo- 
spheric luminosity (L..) and accretion luminosity (Lace), 


L = An Ro TÁ, + GM, Macc /Re- 


In the 2005 epoch, we assume that half of Lao95 = 
1400 Lo arises from a ZAMS photosphere with a so- 
lar metalicity, and the progenitor’s properties are: mass 
M, = 5.4 Mo, radius R, = 2.7 Ro, and effective tem- 
perature Tog = 18000 K (Tout et al. 1996), which lies be- 
tween spectral type B2.5V and B3V (Pecaut & Mamajek 
2013). The accretion rate Mace in the 2005 epoch needed 
to produce the half of Logos is ~ 1.1 x 107° Mo yr7!, 
comparable to the steady accretion rates for the burst- 
ing MYSO of Hunter et al. (2021) and for the 6 Mo 
formation model of Haemmerlé et al. (2013). In the 
2017 epoch, we assume Lace œ~ 8300 Lo, which yields 
Mace © 1.7 x 107? Mo yr-! in the 2017 epoch, where 
M, = 5.4 Mo and R, = 34.1 Ro are assumed. The value 
of Mace in the 2017 epoch estimated from the stellar pa- 
rameters agrees with the value of ~ 2 x 107? Mo yr ^! 
obtained from the proper motions of H20 maser spots 
(Sect. 5.2). In the 2017 epoch, a larger stellar radius 
is required to produce an accretion rate comparable to 
the value estimated from observations. This fact is con- 


sistent with the simulation prediction that a large ac- 
cretion rate may cause the protostellar photosphere to 
bloat, resulting in larger radius (Yorke & Bodenheimer 
2008; Hosokawa & Omukai 2009; Meyer et al. 2019a). 


5.4. Accretion History of M17 MIR 


The amplitude of MIR flux variation is directly pro- 
portional to the enhanced Lace. Nevertheless, the en- 
hanced Lace is a reduced version of enhanced Mace, due 
to the significant inflation of the central object during 
the burst accretion phase. In the burst accretion phase, 
the Mace is amplified by a factor of hundreds compared 
with that in the quiescent phase. M17 MIR even shows 
a WISE/NEOWISE W2 flux variation in a cadence of 
six months, despite the overall trend of becoming con- 
tinuously brighter. The overall trend is suggested to 
originate from the increasing Mace- On the other hand, 
the smaller variation in the cadence of six months might 
imply unsteady Mace With a smaller amplitude. For in- 
stance, the variation in the W2 flux density is —42% be- 
tween 2017.09 and 2018.03. This variation corresponds 
to AMace = —42%, if we assume the stellar radius to 
remain unchanged. This is the largest variation in Macc 
during the burst accretion phase. 

Considering the MIR variation during the decreasing 
phase 1993.03-2004.09, M17 MIR is about three times 
brighter in 1993.03 than in 2017.08. We speculate that 
Mace in 1993.03 is ~ 5 x 107? Mo yr-!, three times that 
in the rebrightening phase. However, Mace in 1993.03 
is probably even higher, since M17 MIR may expand to 
even larger size than that in the rebrightening phase. 

'The MIR light curve of M17 MIR exhibits three phases 
characterized by M17 MIR’s MIR variation in flux den- 
sity over 26 yr. From the above results and analyses, we 
conclude that the decreasing and rebrightening phase of 
M17 MIR indeed reflect the two accretion bursts, which 
are separated by the quiescent phase lasting about 6 yr. 
'The first accretion burst likely began between 1984 and 
1993, and lasted longer than 11 yr but less than 20 yr. 
Thus we use a middle value of 15 yr. We fortunately 
caught the initial stage of the second accretion burst, 
aided by the multiepoch observations of the Spitzer and 
WISE/NEOWISE space telescopes. The duration of the 
second accretion burst is definitely longer than 9 yr. It 
is very likely that M17 MIR has not yet reached to its 
peak brightness, and the duration of the second accre- 
tion burst is similar to that of the first epoch. M17 MIR 
is the first MYSO showing multiple accretion bursts 
with durations of tens of years and a burst magnitude 
of about 2mag. The observed properties of the accre- 
tion bursts of M17 MIR agree with the prediction from 
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Table 9. The Parameters of R17 sp-s-i Models that Satisfy x? — x2,4, < 3Naata for M 17 MIR 
Epoch Model Name x? Av R; Tost Maisk Roe Inclination Comments 
(mag) (Ro) (K) (Mə) (AU) (°) 

(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) 
7Px1sJu9_07 11.387 137.259 8.678 11330.0 0.03119 2707.0 60.37 good fit 
ORya6h0F-07 13.715 138.992 18.77 12710.0 0.01624 1362.0 68.83 L005 &2 L2017 
e6RRFGdF.04 15.246 143.522 8.979 17560.0 0.03541 1359.0 38.46 L2005 © L2017 
1EHMCdxi07 15.652 135.411 17.09 13690.0 0.07064 4221.0 63.84 L2005 © L2017 

2005 e6RRFGdF_05 16.029 143.569 8.979 17560.0 0.03541 1359.0 41.11 L005 © L2017 
wRz4n90OE_04 16.72 120.122 7.19 10480.0 0.02166 3444.0 37.77 good fit 
LwlIekVA6.04 19.687 147.335 16.18 10510.0 0.03299 2554.0 31.59 good fit 
kskLMN4Y_05 20.179 149.157 5.494 24540.0 0.0622 3640.0 47.86 too hot 
kskLMN4Y_06 20.859 148.937 5.494 24540.0 0.0622 3640.0 51.65 too hot 
SUXTjAhLlOS 2.412 134.181 34.14 9631.0 0.07883 4067.0 72.06 good fit 
SUXTjAhl 07 7.531 134.866 34.14 9631.0 0.07883 4067.0 63.34 good fit 

y IvNrh20o.04 10.444 117.76 4.789 22810.0 0.07566 316.5 32.15 bad fit at 19.7, 37.1um 
REwCcjZc.05 10.863 132.046 39.91 8132.0 0.017 4539.0 48.53 bad fit at 19.7, 37.1um 
REwCcjZc.04 11.187 134.554 39.91 8132.0 0.017 4539.0 32.8 bad fit at 19.7, 37.1um 
yIvNrh20o.05 11.204 119.925 4.789 22810.0 0.07566 316.5 42.77 bad fit at 19.7, 37.1um 
2017 REwCcjZc.06 11.619 130.341 39.91 8132.0 0.017 4539.0 56.39 bad fit at 19.7, 37.1um 
REwCcjZc.07 12.488 129.402 39.91 8132.0 0.017 4539.0 60.18 bad fit at 19.7, 37.1um 
V5ttVO0Uz. 07 12.517 145.476 59.43 7264.0 0.0982 4549.0 67.08 bad fit at 19.7, 37.1um 
h2UvYxug.09 13.188 139.651 76.3 9036.0 0.01695 1770.0 82.7 bad fit at 19.7, 37.1um 
LZEBwODZ.08 13.556 138.967 38.38 11690.0 0.01377 941.7 71.79 bad fit at 19.7, 37.1um 
LZEBwODZ.O07 14.313 139.677 38.38 11690.0 0.01377 941.7 60.78 bad fit at 19.7, 37.1um 
jzCvPxfT.06 14.323 126.138 47.92 8185.0 0.0576 4076.0 52.93 bad fit at 19.7, 37.1um 


Table 10. Basic stellar parameters of M17 MIR in the 2005 and 2017 epochs. 


Epoch R. Teg Dx Ay M. Mace Lace 
(Ro) (K) (Lo) (mag) (Me) (Moyr™*) (Lo) 

2005  10.1(3.6) 10863(415) 1400(897) 134.7(10.) 5.4 11x10? 7.0 x 10? 

2017 34.1 9631 9034 134 54 17x10? 8&3x10? 


theoretical simulation that 2mag bursts have a mean 
duration of about 20 yr (Meyer et al. 2019b). 

We assume a bursting accretion timescale of 15 yr, 
which is comparable to the duration of the decreasing 
phase in the MIR light curve (somewhere between 1984 
and 1993 to mid-2004). If M17 MIR can sustain this 
enhanced disk accretion rate of ~ 2 x 107? Mo yr“t, it 
could gain about 0.03 Mo during one accretion burst. 
This amount of accreted mass is an order of magni- 
tude higher than the mass accretion onto S255IR-NIRS3 


(Caratti o Garatti et al. 2017) and two orders higher 
than that onto G358.93-MM1 (Stecklum et al. 2021). 
'The quiescent phase of M17 MIR lasts about 6 yr. Dur- 
ing the time span of 26 yr, the fraction of time in accre- 
tion burst is about 83%. The main accretion phase of 
a high-mass protostellar object in an HMC lasts a few 
times 10* yr (Furuya et al. 2005; Bonfand et al. 2017). 
A crude estimate, assuming tacc = 1 x 10* yr, is that this 
fraction of time in accretion burst yields ~ 14 Mc; mass 
that M17 MIR could accrete from its circumstellar disk. 
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Table 11. Key parameters of six MYSOs with accretion Burst 


Name d LP [burst [99st AFir AF am At Maser Flaring Note 
(kpc) (10°Lo) (10°Lo) (10°Lo) (yr) 
V723 Car 2.5 - 4 ~ 10 - 5 — K 212.9 
S255IR-NIRS3 1.8 29 160 - = 10 2 2 CH30H, H20 K29 
NGC 6334I-MM1 1.3 2.9 47.6 - 16.3 3.90.6 >6 CHa3OH, H2O K-band invisible 
G323.46-0.08 4.8 - — — CHOH K-band visible 
G358.93-MM1 6.75 7.6 19.3 12.7 = 2 e S 0.5 CH30H K 215 
M17 MIR 2.0 1.4 9.0 1.4 >10 = ~2 x15 H20 K > 22 


5.5. Comparisons to other variable MYSOs 


The sample of MYSO accretion bursts is still small. 
Six MYSOs show a luminosity increase due to an ac- 
cretion burst: V723 Car (Tapia et al. 2015), S255IR- 
NIRS3 (Caratti o Garatti et al. 2017; Liu et al. 2018), 
NGC 6334I-MM1 (Hunter et al. 2017, 2021), G323.46- 
0.08 (Proven-Adzri et al. 2019), G358.93-0.03 MM1 
(Stecklum et al. 2021), and M17 MIR in the present 
work. However, these MYSOs are in different pre-main- 
sequence (PMS) evolutionary stages. The durations 
of accretion bursts are poorly constrained except for 
M17 MIR. To qualitatively describe the relation between 
the PMS evolutionary stage and an accretion burst, their 
key parameters are summarized in Table 11. Relatively 
evolved MYSOs are commonly brighter in the near-IR 
K band, while the earliest MYSOs are very faint or 
even invisible in the same band. Generally the K-band 
brightness of these MYSOs is representative of the PMS 
evolutionary stage. Among the six MYSOs tabulated 
in Table 11, four sources show variabilities in the K 
band, indicating relatively late PMS stages. In particu- 
lar, G323.46-0.08 was the first one that showed periodic 
6.7 GHz methanol maser emission with a period of about 
90 days (Proven-Adzri et al. 2019). Later — on the basis 
of multiepoch IR data — it was confirmed as an accretion 
burst MYSO (see Stecklum et al. 2021, and references 
therein). 

For V723 Car, S255IR-NIRS3, and G358.93-MM1 
the duration of the accretion burst At was constrained 
from their IR luminosity changes. The durations of 
their luminosity increases are of order of years and thus 
comparable to each other. The two earliest MYSOs, 
NGC 6334I-MM1 and M17 MIR, are both embedded in 
HMCs. The progenitor of NGC 6334I-MM1 in its pre- 
burst phase might be able to produce a hypercompact 
region MMIB (Brogan et al. 2016), implying that NGC 
6334I-MM1 is more massive and at a later evolutionary 
stage than M17 MIR. Acoording to the schematic evo- 
lutionary sequence of a massive protostar (Motte et al. 
2018), M17 MIR might be the earliest MYSO with ac- 
cretion outbursts. Considering the extreme youth of 


M17 MIR, its multiple accretion bursts suggest that mi- 
nor accretion bursts are frequent in the very early stages 
of massive star formation. 


6. CONCLUSIONS 


'The multiwavelength and multiepoch data of the mas- 
sive protostar M17 MIR and the subsequent analyses of 
these data yield the following major results. 

1. By analyzing the IR data of M17 SW in the 
3 — 70 um range, M17 MIR is identified as the driv- 
ing source of the expanding structure traced by H20 
maser motions. The sub-mm dust continuum emission 
and molecular line emission indicates that M17 MIR is 
embedded within a HMC. 

2. The unified flux density of M17 MIR at the 3.8 um 
L’ band and the 4.6 um WISE W2 band clearly shows 
three stages in the light curve: i) the decreasing phase 
during 1993.03 to mid 2004 with an amplitude of 10 
at 3.8 um, ii) the quiescent phase during mid 2004 to 
mid 2010, and iii) the rebrightening phase since mid 
2010 until now with an amplitude of 5 at 4.6 jum. The 
22 GHz H2O maser emission associated with M17 MIR 
also shows flux variations. They are correlated with the 
MIR flux variation in time space, which together indi- 
cate enhanced disk accretion rate onto M17 MIR during 
the decreasing and rebrightening phase. 

3. The kinematics of H2O maser spots, measured in 
2009.10—2010.12, in the close vicinity of M17 MIR yield 
a disk accretion rate Mace ~ 7 x 1074 Mo yr-!. In the 
later stage of the rebrightening phase after 2014, Macc is 
~ 2 x 107? Mo yr- 1, given by the 4.6 um flux increased 
by a factor of ~ 3 compared to the initial stage in 2010. 

4. Radiative transfer modeling for the SEDs in the 
2005 and 2017 epoch constrains the basic stellar param- 
eters of M17 MIR. In the quiescent phase in the 2005 
epoch, M17 MIR likely has total luminosity Loo95 = 
1400 + 897 Lo, stellar mass M, ~ 5.4 Mo, stellar ra- 
dius R, = 10.1 + 3.6 Rọ, and accretion rate Macc ~ 
1.1 x 107? Mo yr7!. In the accretion burst phase (e.g. 
in the 2017 epoch), L217 = 9034 Lo5( AL ~ 7600 Lo), 
R, = 34.1 Ro, and Mace ~ 1.7 x 107? Mo yr7!. In the 
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accretion outburst, a larger stellar radius is required to 
produce Macc consistent with the value estimated from 
the kinematics of H20 maser spots. 

5. The decreasing and rebrightening phase of 
M17 MIR reflect two accretion bursts (At ~ 9 — 20 yr), 
which are separated by the quiescent phase lasting about 
6yr. During the time span of 26 yr, the fraction time 
of accretion burst is about 88%. M17 MIR is the first 
discovered massive protostar showing multiple accretion 
bursts with durations of tens of years and a burst mag- 
nitude of about 2mag. Among the limited sample of 
accretion burst MYSOs, the extreme youth of M17 MIR 
suggest that minor accretion bursts are frequent at very 
early stages of massive star formation. 

6. Even during the current accretion burst (since mid 
2010 to until now), the W2 flux density of M17 MIR 
shows variations of small amplitude up to a factor of 3, 
indicating variable accretion rate with amplitutde much 
lower than the overall contrast between queiscent and 
outburst state. Long time monitoring at MIR bands will 
record the bursting accretion activity of M17 MIR in the 
near future. Because M17 MIR is still in its earliest stage 
of evolution, it will be a unique testbed for probing the 
earliest phases of massive star formation. 
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